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Aerodynamical ly  stable  operation  of  the  compression  system 
is  vital  in  an  aircraft  turbine  engine.  A  compressor  delivers  a 
desired  quantity  of  airflow  at  a  desired  pressure  level  during 
the  stable  mode  of  operation.  The  aerodynamic  stability  limit 
line  or  surge  line  is  the  critical  line  on  a  compressor  map  which 
is  a  plot  of  the  total  pressure  ratio  across  the  compressor  versus 
corrected  airflow  rate.  Only  below  the  surge  line  is  the  com¬ 
pressor  stable.  Compressor  surge  adversely  affects  engine  per¬ 
formance  by  reducing  engine  airflow  and  thrust.  The  surge  line 
is  usually  determined  experimentally,  based  on  steady-state  and 
some  limited  unsteady  data  and  various  empirical  correlations 
and  corrections.  The  surge  line,  however,  is  affected  by  the 
fluctuations  and  distortions  of  the  inlet  flow  and  other  transient 
phenomena  during  the  operation  of  the  system.  Since  it  is  not 
feasible  to  experimentally  evaluate  a  compressor  or  an  engine 
under  all  possible  conditions,  mathematical  models  are  built  to 
simulate  the  operation  of  compression  systems. 

Kimzey^ developed  a  one -dimensional  time -dependent  model 
for  the  analysis  of  the  effects  of  planar  disturbances  on  a 
compressor  system  on  the  basis  of  conservation  laws  of  mass, 
momentum  and  energy.  Experimental  data  are  used  to  synthesize 
the  stage  characteristics  of  the  compressor.  The  compressor  stage 
force  and  shaft  work  which  are  needed  in  the  model  are  calculated 
based  on  these  characteristics.  Other  factors  such  as  the  force 
of  compressor  casing  acting  on  the  fluid,  heat  added  to  the  fluid 
and  compressor  bleed  flow  rate  are  also  included  in  the  model 
equations.  The  conservation  equations  are  discretized  spatially 
by  the  use  of  a  two-sided  difference  scheme.  Boundary  conditions 
are  imposed  on  total  pressure  and  total  temperature  at  the  inlet 
boundary  and  on  static  pressure  or  airflow  rate  at  the  exit  bound¬ 
ary.  The  resulting  system  of  ordinary  differential  equations  are 
integrated  forward  in  time  by  a  fourth-order  Runge-Kutta  scheme. 
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These  models  have  been  applied  to  a  variety  of  compression 
systems  by  Kimzey^^  and  Chamblee,  Davis  and  Kimzey^),  They 
are  used  to  analyze  and  extrapolate  the  experimental  data  and  to 
study  the  effects  of  unsteady  disturbances  on  the  aerodynamic 
stability  of  the  compression  system.  These  models  are  not  always 
numerically  stable.  Some  of  the  techniques  used  for  overcoming 
the  numerical  instabilities  are  to  increase  the  friction  coeffi¬ 
cient  in  the  inlet  and  exit  ducting,  altering  the  ducting  lengths 
or  areas  and  time  averaging  the  numerical  solutions  once  every 
few  time  steps.  However,  these  techniques  can  achieve  only 
conditional  numerical  stability  for  some  of  the  compression 
systems.  Also  the  application  of  the  models  has  been  restricted 
to  limited  regions  of  the  operating  map.  For  certain  regions  of 

the  operating  map  the  models  are  numerically  unstable. 

(2) 

Davis  v  '  has  applied  the  MacCormack's  explicit  finite 
difference  scheme  to  solve  the  partial  differential  equations  of 
the  model  and  an  approximate  version  of  the  method  of  character¬ 
istics  to  impose  the  boundary  conditions.  This  scheme  is  more 
stable  numerically  than  the  earlier  method,  but  this  also  exhibits 
numerical  oscillations  in  the  solutions  under  certain  conditions. 
These  oscillations  are  controlled  or  reduced  by  the  addition  of 
extra  friction  or  dissipation  in  the  inlet  duct.  However,  the 
additional  extra  friction  may  degrade  the  accuracy  of  the  simu¬ 
lation  of  the  actual  physical  system. 

In  this  report  we  present  a  study  of  the  numerical  stability 
of  the  compressor  model  and  develop  different  schemes  of  enhanced 
numerical  stability.  The  following  items  are  discussed. 

a.  mathematical  study  of  the  compressor  model  equations 
and  the  boundary  conditions  (application  of  the 
method  of  characteristics  for  imposing  the  boundary 
conditions  accurately) 

b.  numerical  schemes  and  their  stability  characteristics: 

1.  Runge-Kutta  Scheme 

2.  JRS  (Jameson,  Rizzi,  Schmidt)  Scheme 

3.  MacCormack  Scheme 
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c.  stable  test  calculations  using  the  AEDC  data 
The  alternate  methods  of  solution  for  the  compressor  model 
equations  developed  in  this  report  have  been  shown  to  be 
numerically  stable  and  are  compared  against  each  other  and  with 
some  limited  experimental  data. 
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II.  COMPRESSOR  MODEL  AND  THE  BOUNDARY  CONDITIONS 


The  compression  system  modeled  and  the  control  volume 

system  are  shown  in  Fig.  1.  The  governing  equations 
q  2) 

derived  v  '  by  the  application  of  the  mass,  momentum  and 
energy  conservation  laws  to  the  elemental  control  volume  in 
Fig.  1-d  in  which  blade  forces  ,  wall  shear  forces,  shaft 
work  done,  heat  added  to  the  fluid  and  mass  bleed  flow  rate  are 
included.  The  resulting  system  of  first  order  partial  differen¬ 
tial  equations  can  be  written  as  follows: 


+  ||(3)  +  g(u,z,t)  = 


9z 


(1) 


where 


pA 

pAU 

->  ■+ 

9 

pAU 

,  f  (u)  = 

A  (pU  +  P) 

2 

_PA  (e  + 

- PAUCpTt  - 

g(u,5?,t)  = 


W, 


B 


_p 

9z 


-W 


(2) 


The  various  symbols  represent  the  following: 
p  =  density 
A  =  area 

U  =  axial  velocity 
e  =  internal  energy 
P  =  static  pressure 

C  =  specific  heat  at  constant  pressure 
P 

T^_  =  stagnation  temperature 
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=  compressor  bleed  flow  rate 
F  =  force  of  compressor  blading  and  casing  friction 
acting  on  the  fluid 

Wg  =  stage  shaft  work  added  to  fluid  in  control  volume 
Q  =  rate  of  heat  addition  to  control  volume 

In  conjunction  with  the  partial  differential  equations  (1) , 
we  have  the  equation  of  state  given  by 

P  =  p  R  T  (3) 

where  T  is  the  static  temperature  and  R  the  gas  constant.  In 
addition  we  have  the  equations  relating  the  internal  energy  and 
stagnation  temperature  to  other  variables  as  follows: 


e  =  C  T  /y 

p  ’ 

(4) 

C  T^  =  C  T  +  — 

p  l  p  2 

(5) 

The  area  distribution  A  (z)  is  known  for  a  given  compressor 
system.  In  the  vector  g,  which  acts  like  a  forcing  function  in 
the  differential  equation,  various  terms  are  modeled  empirically 
for  a  particular  compressor  except  perhaps  the  pressure  term 
which  may  be  computed  as  a  part  of  the  solution.  F(z,t)  represents 
the  forces  of  the  blades  and  the  casing  friction  acting  on  the 
fluid.  In  practice  it  is  difficult  to  isolate  F(z,t)  empirically 
from  the  experimental  data  and  hence  the  total  term  (F  +  P  8A/3z) 
which  represents  the  forces  including  the  wall  pressure  area 
force  is  modeled  from  the  experimental  data.  (z,t)  is  the 

O 

shaft  work  done  on  the  fluid,  which  depends  on  the  stage  character¬ 
istics  of  the  compressor.  The  stage  characteristics  are  modeled 
empirically  based  on  experimental  measurements  of  stage  total 
temperature,  flow  rate,  total  pressure  and  flow  angularity  at  the 

stage  entry  and  exit,  and  are  corrected  to  account  for  unsteady 

(4) 

flows  through  cascades  based  on  the  work  of  Goethert  and  Reddy  . 
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Q(z,t)  is  tiie  heat  addition  to  the  fluid  and  W^Cz.t)  is  the 

D 

bleed  flow  distribution  function  which  models  the  mass  removal 
or  addition  between  stages. 

The  equations  (1)  can  be  written  in  the  quasi-linear  form 

If  +  M(3)  ll  +  g  (u,z,t)  =  0  (6) 

where  M(u)  =  —  is  called  the  Jacobian  Matrix  of  the  system.  As 
3u 

it  is  shown  in  the  Appendix  A,  the  Jacobian  matrix  M  has  three 
real  and  distinct  eigenvalues  U,  U+c ,  U-c,  where  c  is  the  speed 
of  sound.  Hence  the  system  of  partial  differential  equations 
(1)  is  hyperbolic.  According  to  the  theory  of  hyperbolic  partial 
differential  equations ,  the  number  of  boundary  conditions  to 
be  imposed  on  the  left  boundary  should  equal  the  number  of  positive 
eigenvalues  of  M,  which  in  the  present  case  is  two  for  subsonic 
incoming  flows.  In  other  words  there  are  two  incoming  character¬ 
istics  at  the  inlet  for  subsonic  flows  and  correspondingly  we 
need  to  impose  two  boundary  conditions  there.  Similarly  on  the 
right  boundary,  the  number  of  boundary  conditions  to  be  imposed 
should  equal  the  number  of  negative  eigenvalues  of  M  there.  For 
subsonic  outflows,  as  is  the  case  in  the  present  model,  we  have 
one  negative  eigenvalue,  U-c,  which  corresponds  to  one  incoming 
characteristic  at  the  outflow  boundary  and  we  impose  one  right 
boundary  condition.  In  the  present  model,  we  prescribe  the  total 
pressure  Pt  (t)  and  the  total  temperature  Tt  (t)  at  the  left 
boundary  and  the  static  pressure  P  (t)  or  the  mass  flow  rate  W  (t) 
at  the  right  boundary. 

Since  it  is  a  hyperbolic  system,  the  best  way  to  formulate 
three  other  conditions  at  the  boundaries,  which  may  be  necessary 
in  numerical  computations  is  by  the  method  of  characteristics. 

The  following  three  compatibility  equations  along  the  three 
characteristic  directions  are  derived  in  Appendix  A.  These 
equations  are  equivalent  to  the  original  system  of  partial 
differential  equations. 
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„2  dp  dP  ,  2  ~  ~  p.  -t  dz  TT 

c  at  -  at  +  c  gi  _  g3  =  0  alons  at  =  u  (7) 

pc  ai  +  af  +  +  §3  =  0  alons  af  =  u+c  (8) 

-P°  +  g|  -  P°  g2  +  §3  =  0  along  =  U-c  (9) 

where  the  quantities  g^  ,  g£ ,  §3  are  defined  in  the  Appendix  A. 

The  first  two  characteristic  curves,  corresponding  to  the 
eigenvalues  U  and  U+c  are  the  so-called  right-running  character¬ 
istics  and  the  third  one  corresponding  to  U-c  is  the  left-running 
characteristic.  We  have  two  boundary  conditions  and  one  compati¬ 
bility  equation  (9)  at  the  inflow  boundary  and  one  boundary  condi¬ 
tion  and  two  compatibility  equations  (7)  and  (8)  at  the  outflow 
boundary.  These  equations  are  sufficient  to  solve  for  all  the 
flow  variables  at  the  boundaries  necessary  for  the  interior  point 
numerical  schemes  considered  in  this  report.  An  accurate  itera¬ 
tive  method  is  used  to  solve  these  equations  on  the  boundaries. 
Details  are  given  in  the  Appendix  B. 


13 


AEDC-TR-82-1  6 


III.  NUMERICAL  SCHEMES 


Three  different  numerical  schemes,  namely  the  Runge-Kutta 
Scheme,  the  JRS  (Jameson,  Rizzi  and  Schmidt)  Scheme  and  the 
MacCormack  Scheme  have  been  used  to  integrate  the  partial 
differential  equations  forward  in  time.  All  of  the  schemes  are 
explicit.  In  the  JRS  and  the  Runge-Kutta  Schemes,  the  spatial 
derivatives  are  approximated  by  second-order  accurate  two-sided 
differences,  while  in  the  MacCormack  Scheme  forward  and  back¬ 
ward  differences  are  used  in  two  steps  which  also  yield  second 
order  accuracy. 

A.  RUNGE-KUTTA  SCHEME 

The  partial  differential  equations  (1)  can  be  rewritten  as 


(-)  • 
^  9t'  i 


=  -c  + 

^  9z  ^ 


g) 


(10) 


where  i  is  the  interior  point  index.  The  spatial  derivatives 
on  the  right  are  approximated  by  second-order  accurate  two-sided 
differences .  The  resulting  set  of  ordinary  differential  equations 
are  of  the  form: 


du  |  _ 

cTt '  i 


t) 


i  =  2, . . .1-1 


(11) 


where  the  index  i  =  1  corresponds  to  the  left  boundary  and  the 
index  i  =  I  to  the  right  boundary. 

The  equations  (11)  are  integrated  forward  in  time  by  a 
fourth  order  accurate  (in  time)  Runge-Kutta  Scheme: 


+n+l 

u 


±n  ,  1 

u  +  T 


(k-i 


+  2tr 


O'v 


where  1-^  =  At  F  (un,tn) 

^2  =  At  F  (un  +  lc^/2,  tn  +  At/2) 
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L  =  At  F  (un  +  &0/2,  t  +  At/2) 
j  n 

^  =  At  F  (un  +  ^3 ,  tn  +  At) 


The  local  linearized  stability  analysis  (von-Neumann 
Stability)  for  the  difference  scheme  (11)  and  (12)  is  carried 
out  as  follows.  Since  the  non-homogeneous  terms  g  in  the 
equations  (1)  or  (6)  are  computed  at  each  time  step  using  some 
empirical  models  based  on  experimental  data,  we  need  only 
consider  the  homogeneous  part 


9u 

at 


+  M 


3u 

az 


0 


(13) 


where  M  =  9f/9u  is  the  Jacobian  matrix,  for  the  purpose  of 
numerical  stability  analysis.  Even  though  we  use  non-uniform 
spatial  steps  in  the  numerical  computations,  we  consider  the 
central  difference  scheme  with  uniform  mesh  steps  in  the  z-direc- 
tion  and  constant  (locally  frozen)  matrix  M  for  the  purpose  of 
analysis.  We  can  approximate  (13)  by 


_  „  *1+1  A-l 

at  u  2AX 


(14) 


-  "  Iky 

We  analyze  how  a  Fourier  wave  component  of  the  form  A^(t)e  “  , 
where  I  =  -yJ-T  and  k  is  the  wave  number,  is  propagated  in  time 
by  the  difference  equations.  For  a  stable  numerical  scheme  the 
magnitudes  of  the  wave  amplitudes  should  remain  bounded  for  all 
frequencies  of  interest.  Otherwise  the  numerical  scheme  is 
considered  unstable.  The  time  growth  of  a  Fourier  component  of  the 
solution  with  wave  number  k  is  obtained  by  seeking  the  solution 
of  the  equation  (14)  in  the  form 


3±  =  K  (t)  eIklAx 


v.  (_ ;  c  ,  where  X  —  1 

Substituting  u^  in  the  equation  (14)  we  obtain 


(15) 


dt 


I  M  SinlkAxW 
Ax  k 


(16) 
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where  o  <  k  Ax  <  it.  The  limitation  k  Ax  <  tt  is  imposed  by  the 
fact  that  lc  =  tt/ A x  is  the  highest  wave  numbei:  that  can  be  resolved 
by  a  finite  difference  scheme  with  mesh  width  Ax.  The  time  growth 
of  the  amplitudes  is  thus  governed  by  a  linear  system  of 
ordinary  differential  equations  (16).  The  coefficient  matrix, 

B  =1  MSin(kAx)/Ax  of  this  system  has  three  eigenvalues  equal  to 
I  X  Sin(kAx)/Ax,  where  X  represents  the  three  eigenvalues  of  M, 
namely,  U,  U-c  and  U+c.  Since  the  eigenvalues  of  the  coefficient 
matrix  lie  on  the  imaginary  axis  of  the  complex  plane,  the  system 
of  ordinary  differential  equations  (16)  is  neutrally  stable. 

When  a  Runge-Kutta  Scheme  is  used  to  integrate  the  system 
of  ordinary  differential  equations  (16)  ,  we  can  obtain  the  stabil¬ 
ity  condition  from  the  theory  of  numerical  integration  of  ordinary 
differential  equations . The  eigenvalues,  y  of  the  matrix  At  B 
(At  times  the  coefficient  matrix)  must  lie  within  a  zone  of 
stability  in  the  complex  plane  corresponding  to  the  particular 
Runge-Kutta  Scheme.  The  eigenvalues  y  are  given  by 

y  =  I  X  ^  Sin(kAx)  ,  o  <  kAx  <  tt  (17) 

where  X  are  the  eigenvalues  of  M,  namely,  U,  U-c,  U+c.  The 
stability  zone  in  the  complex  plane  corresponding  to  the  fourth 
order  Runge-Kutta  Scheme  is  shown  in  Figure  2.  The  upper  most 
point  on  the  imaginary  axis  within  the  region  of  stability  is 
approximately  2.7.  Since  all  the  three  values  of  y  lie  on  the 
imaginary  axis,  the  condition  for  stability  is 

Max  |  X  ^  Sin  kAx  |  <2.7  ,  o  <  kAx  <  it 

This  gives  the  necessary  condition  for  stability  in  terms  of  the 
Courant  number,  CN  as 


CN  =  Max  (U+c)  ^  <2.1  (18) 

Even  if  the  time  step  At  is  chosen  to  satisfy  the  equation  (18) , 
stability  is  still  marginal  because  the  non-uniformness  of  the 
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mesh  and  the  non-linearities  in  the  convection  terms  can  make  y 
to  lie  in  the  unstable  region  to  the  right  of  the  imaginary  axis 
of  the  complex  plane.  In  fact,  non-linear  instabilities  do  occur 
with  the  Runge-Kutta  scheme  in  practice.  In  the  actual  numerical 
calculations  it  has  been  found  necessary  to  time  average  the 
solution  once  every  5  or  10  time  steps  to  suppress  the  non-linear 
instabilities  even  after  obeying  the  Courant  number  restriction. 
Those  details  are  discussed  in  the  section  dealing  with  the 
results . 

B.  JRS  (JAMESON,  RIZZI  AND  SCHMIDT)  SCHEME 

This  scheme  proposed  by  Jameson,  Rizzi  and  Schmidt^^  is 
also  based  on  two-sided  spatial  difference  approximations  of 
second  order  accuracy  similar  to  the  Runge-Kutta  Scheme.  It  is 
a  three-step  method  with  one  predictor  and  two  corrector  steps. 

Let  D  denote  a  second  order  accurate  two-sided  difference 
approximation  for  9/3z.  The  JRS  Scheme  can  be  written  as  follows: 


P: 

->n 
=  ui 

Cl: 

3J+1 

+n 
=  u . 

l 

C2: 

+n+l 

ui 

->-n 
=  ui 

-  At  (D*“  +  IJ) 

-  i  At  ml  +  1"  +  Dlfr  +  if1) 

-  —  At  (dJ?  +  +  DiV+1  +  g?+1) 

2  v  l  6i  l  &i  ' 


(19) 


(7) 

The  local  stability  analysis v  '  (von-Neumann  stability 
analysis)  of  this  scheme  can  be  carried  out  by  analyzing  the 
growth  of  the  amplitudes  of  the  basic  Fourier  modes  of  the 
numerical  solution  from  one  time  step  to  the  next.  Since  the 
discretization  in  space  and  time  are  combined  together  in  equation 
(19) ,  we  seek  solutions  of  the  form 


=  tn  PIkiAx 

Ui  k  e 


,  where  I  = 


(20) 
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For  the  purpose  of  linear  stability  analysis,  the  non-linear  term 
Df^1  is  replaced  by  M  Du^1  where  M  is  the  Jacobian  matrix  (frozen 
locally)  and  D  is  a  central  difference  operator.  Also  the 
nonhomogeneous  terms  g  are  dropped  as  mentioned  before,  since 
they  do  not  influence  the  stability  analysis.  After  some  algebra, 
one  can  obtain  a  relation  between  the  amplitudes  at  the  nth 
and  (n+1)  th  time  steps  in  the  form 

££+1  =  g  ££  (21) 

where  G  is  known  as  the  amplification  matrix.  It  can  be  shown 
that 


G  =  (I  -  B  +  y  B2  -  B3)  (22) 

where  I  is  the  identity  matrix  and  B  =  I  At  M  Sin(kAx)/Ax.  Linear 
stability  is  assured  if  Max | eigenvalues  of  G|<1  . 

Eigenvalues  of  G  =  1  -  y  +  —  (23) 

where  y  is  given  by  the  equation  (17) . 

o  -i  /  1 

|  Eigenvalues  of  G | “  =  1  -  v  (1  -  v  ) 

<1  ,  if  | v  |  <2 

where  v  =  A  At  Sin(kAx)/Ax  and  X  =  U,  U-c  or  U+c . 
condition  for  linear  stability  is  obtained  as 

CN  =  Max  (U+c )  <2  (25) 

The  time  accuracy  of  the  scheme  is  of  second  order,  which  is  less 
than  that  of  the  Runge-Kutta  Scheme.  However  it  requires  less 
computing  per  time  step  while  the  time  step  can  be  taken  nearly 
as  large  as  that  of  the  Runge-Kutta  Scheme.  However,  this  is  also 
a  non-dissipative  scheme  just  like  the  Runge-Kutta  Scheme  and  non¬ 
linear  instabilities  occur  under  certain  conditions.  For  certain 
test  runs  it  has  been  found  necessary  to  time  average  the  solutions 
once  every  5  or  10  time  steps . 


(24) 

Thus  the 
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C.  MacCORMACK  SCHEME 

One  of  the  most  popular  schemes  for  solving  the  time  depen¬ 
dent  Euler  equations  or  Navier-Stokes  equations  is  the  MacCormack 
Scheme^)  due  to  its  many  desirable  features.  It  is  a  two-step, 
predictor-corrector  method.  In  the  predictor,  the  spatial 
derivatives  are  approximated  by  a  backward  difference  scheme  and 
in  the  corrector  by  a  forward  difference  scheme.  Although  either 
the  predictor  or  the  corrector  is  only  first  order  accurate  in 
time  and  space,  the  total  scheme,  however,  is  second  order 
accurate  both  in  time  and  space.  The  scheme  can  be  written  as 


+n+T 
u . 

1 


+n+l 

ui 


+n 
u . 
1 


_1_ 

2 


IT.  -  h-d  -  At  Si 


+n  ,  ->n+l 
ui  +  ui 


At  ^jn+1  _  jn+r} 
Azi+1  1+1  i 


(26) 

_ \ 

.  .  ->n+l 

At  gi 


where  Az^  =  zp“zp„q-  Using  an  analysis  similar  to  that  outlined 
for  the  JRS  scheme  one  can  show  that  the  necessary  condition  for 
local  (linearized)  stability  of  this  scheme  is  that  the  Courant 
number,  CN  obey 


CN  =  max  (U+c)  <  1  (27) 

MacCormack  scheme  is  dissipative  in  the  sense  that  the  high 
frequency  components  of  the  errors  are  dissipated  rapidly.  Thus 
it  is  less  susceptible  to  non-linear  instabilities  than  the  Runge- 
Kutta  or  JRS  schemes.  As  a  matter  of  fact  whenever  the  Courant 
number  is  less  than  1,  the  numerical  calculations  have  always 
been  stable  without  any  time  averaging.  While  the  time  step 
that  can  be  chosen  for  this  scheme  is  less  than  those  of  the 
Runge-Kutta  and  JRS  schemes,  the  computation  per  step  is  less. 

With  the  three  difference  schemes  discussed  above,  the 
boundary  conditions  are  imposed  using  the  method  of  characteris¬ 
tics  as  discussed  in  the  last  section  and  the  Appendix  B.  Test 
calculations  have  been  performed  using  all  the  schemes  and  the 
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results  are  discussed  in  the  next  section.  Simple  extra¬ 
polation  techniques  at  the  boundaries  have  been  found  to  be 
either  unstable  or  otherwise  unsatisfactory. 
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IV.  TEST  CALCULATIONS 

y 

The  model  described  in  the  last  two  sections  has  been 

applied  to  a  typical  high-pressure  compression  system.  A 

Schematic  of  this  system  is  presented  in  Fig.  1-c  which 

includes  4  inlet  ducting  control  volumes,  10  compressor  stage 

control  volumes  and  5  combustor  control  volumes.  For  the  inlet 

ducting  and  combustor  control  volumes  the  friction  force  is 

2 

calculated  by  =  cf  pv  / 2  where  cf  is  the  assumed  skin  fric¬ 
tion  factor.  The  compressor  blade  force  and  the  shaft  work  in 
each  compressor  stage  are  calculated  from  the  corresponding 
stage  characteristics v  '  .  Compressor  bleed  flow  rate  and  heat 
addition  to  the  fluid  are  not  considered  in  the  test  cases  run. 

In  the  boundary  conditions,  the  total  pressure  and  the  total 
temperature  at  the  inlet  are  fixed,  but  the  static  pressure  at 
the  exit  is  specified  as  a  function  of  time.  The  exit  pressure 
is  increased  with  time  as  a  linear  function  (pressure  ramp)  at 
50  psia/sec  in  one  case  and  20  psia/see  in  another  case.  Com¬ 
putations  have  shown  that  the  compressor  surge  occurs  at  about 
the  same  exit  pressure  level,  as  shown  in  Fig.  3.  Hence^  most 
of  the  computations  have  been  done  with  50  psia/sec  exit  pressure 
ramp . 

For  each  run,  the  initial  values  of  the  dependent  variables 
namely  the  density,  the  mass  flow  rate  and  the  total  energy  are 
input  at  all  the  spatial  locations  of  the  grid.  They  can  be 
obtained  from  another  steady  state  code  or  can  be  empirically 
guessed.  Force  calculation  data  such  as  the  stage  characteristics 
are  also  input  for  a  constant  corrected  speed.  Time  dependent 
calculations  are  done  for  a  fraction  of  a  second  (0.1  sec  in  the 
present  case)  with  fixed  boundary  conditions  at  both  ends  in 
order  to  reach  a  steady-state  starting  solution.  At  this  point 
the  exit  static  pressure  is  increased  at  a  rate  of  50  psia/sec. 
Test  runs  have  been  made  with  two  different  corrected  speeds; 
one  at  87  percent  of  design  corrected  speed  and  the  other  at  102 
percent  corrected  speed. 
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A  test  case  with  87  percent  corrected  speed  has  been  run 
with  At  =  0.00005  which  corresponds  to  the  Courant  number, 

CN  —  0.6  using  all  three  difference  schemes.  The  Runge-Kutta 
scheme  requires  time -aver aging  of  the  solution  once  every  5  time 
steps  in  order  to  suppress  the  non-linear  numerical  instabilities. 
The  numerical  solution  obtained  by  the  JRS  scheme  exhibits  some 
numerical  oscillations  but  they  are  not  severe  enough  to  cause 
non-linear  instability.  However  these  numerical  oscillations 
are  removed  by  time-averaging  the  solution  once  every  5  time  steps 
in  the  JRS  scheme.  The  MacCormack  Scheme  produces  stable  solution 
without  any  numerical  oscillations  till  the  point  of  surge  and 
does  not  require  any  time -averaging .  The  pressure  ratio  across 
the  compressor  is  plotted  against  the  percent  corrected  compressor 
inlet  airflow  in  Figs.  4,  5,  and  6  corresponding  to  the  Runge- 
Kutta,  JRS  and  MacCormack  Schemes  respectively.  All  these  schemes 
produce  the  same  performance  map  except  near  the  surge  point.  The 
Runge-Kutta  and  JRS  Schemes  predict  the  surge  point  within  2  per¬ 
cent  of  the  pressure  ratio  of  the  experimental  surge  point.  The 
MacCormack  scheme  does  it  even  better  and  predicts  the  surge  point 
within  1  percent.  Figures  7,  8  and  9  show  the  pressure  ratio 
across  the  compressor  and  the  percent  corrected  airflow  against 
time  till  the  point  of  surge, obtained  by  the  three  methods  and 
they  compare  very  well  against  each  other. 

In  addition  to  the  overall  performance  plots,  we  have  also 
plotted  some  computed  results  of  the  individual  stages.  Figures 
10,  11  and  12  show  individual  stage  entrance  static  pressures 
against  time  and  Figs.  13,  14  and  15  show  the  stage  entrance 
airflow  rate  against  time  calculated  by  the  three  methods,  just 
before  the  flow  breaks  down.  At  the  time  of  compressor  surge, 
the  stage  that  initiates  the  surge  will  indicate  a  sudden  increase 
in  stage  entrance  pressure  and  a  sudden  decrease  in  exit  pressure 
(which  is  the  entrance  pressure  of  the  next  stage) .  This  signa¬ 
ture  is  caused  by  flow  stagnation  at  the  entrance  and  recirculation 
at  the  stage  exit,  perhaps  caused  by  flow  separation  on  the  blades. 
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By  examining  the  individual  stage  pressure  performance  in  the 
Figs.  10,  11  and  12,  it  may  be  inferred  that  the  2nd  and  7th 
stages  initiate  the  surge  process.  Experimental  data,  however, 
indicates  that  the  surge  is  initiated  at  the  7th  stage.  It 
seems  that  the  stage  characteristic  data  for  the  2nd  stage,  used 
in  this  test  run,  is  not  very  accurate. 

The  model  predictions  of  the  compressor  performance  for  the 
case  of  102  percent  corrected  speed  are  given  in  Figs.  16  to  27. 

The  experimental  surge  point  for  this  case  is  at  100.4  percent 
corrected  compressor  inlet  airflow  when  the  pressure  ratio  across 
the  compressor  equals  to  10.  All  three  methods  predict  the  surge 
point  within  0.5  percent  of  the  experimental  values.  Figures  22 
to  24  indicate  that  the  surge  process  is  initiated  at  the  7th 
stage  in  this  case.  In  general  the  model  predictions  in  this 
case  by  the  three  methods  agree  with  each  other  well  and  are 
better  than  the  previous  case.  Experimental  data  agrees  with  the 
prediction  of  the  7th  stage  initiating  the  surge  in  this  case. 

We  have  made  computations  with  larger  time  step  (At  =  0.0001) 
which  corresponds  to  a  Courant  number  CN  —1.2  to  see  if  the  runs 
can  be  speeded  up.  The  Runge-Kutta  and  JRS  Schemes  with  time- 
averaging  yield  same  results  as  before  for  both  87  percent  and 
102  percent  cases,  but  require  only  half  the  computing  time.  The 
MacCormack  Scheme  which  is  stable  only  for  CN  <  1 ,  works  even  with 
CN  —  1.2  for  the  87  percent  case  and  gives  the  same  results  as  befor 
in  half  the  computing  time.  The  MacCormack  Scheme  with  CN  —  1.2 
becomes  unstable  for  the  102  percent  case,  but  can  be  stabilized 
however,  by  time-averaging  the  solution  once  every  2  time  steps 
in  which  case  it  yields  the  same  results  as  before  again  in  half 
the  computing  time.  In  general  it  may  be  better  to  use  the 
MacCormack  method  with  CN  —  1.0  and  no  time-averaging. 

In  conclusion  the  Runge-Kutta  and  JRS  Schemes  with  time¬ 
averaging  and  the  MacCormack  Scheme  without  any  time -aver aging 
and  the  boundary  conditions  imposed  with  the  use  of  characteris¬ 
tics  in  all  the  schemes  yield  numerically  stable  solutions  of  the 
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compressor  model.  Among  the  three  schemes  the  MacCormack  Scheme 
is  preferred  for  its  stability  characteristics  and  efficiency  in 
computation. 
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V.  CONCLUSIONS 

A.  The  one-dimensional  compressor  model  equations  are  hyper¬ 
bolic  and  the  imposition  of  proper  boundary  conditions 

is  critical  for  solving  them.  For  obtaining  the  numerical 
solutions  of  these  equations,  imposition  of  the  boundary 
conditions  based  on  the  method  of  characteristics  has  been 
found  to  be  the  most  accurate  and  stable  technique. 

B.  Three  different  finite  difference  schemes  have  been  analyzed. 
The  Runge-Kutta  Scheme  and  the  JRS  (Jameson,  Rizzi  and 
Schmidt)  Scheme  exhibit  non-linear  instabilities  even  when 
the  Courant  number  restriction  is  obeyed.  However,  both 

the  schemes  can  be  stabilized  by  time-averaging  the  solu¬ 
tions  every  5  time  steps.  MacCormack  Scheme  has  been 
found  to  yield  numerically  stable  solutions  without  any 
t ime -aver ag ing . 

C.  All  of  the  three  difference  schemes  together  with  the 
boundary  conditions  based  on  the  method  of  characteristics 
predict  compressor  surge  reasonably  accurately  and  agree 
with  experimental  data  reasonably  well  for  the  test  cases 
run.  Overall,  MacCormack  Scheme  has  been  found  to  be  more 
accurate  and  reliable  than  the  other  two  schemes. 
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APPENDIX  A 

DERIVATION  OF  THE  CHARACTERISTICS 
AND  COMPATIBILITY  EQUATIONS 


The  one -dimensional  time-dependent  compressor  model  equations 
in  the  conservation  form,  as  given  in  equation  (1)  are 


9u(z ,  t) 
9 1 
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3f 

9z 


(u) 


+  g  (u ,  z  ,  t )  =  0 
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pAU  (e  +  y  + 

_  pU(h  +  ) 

g(u,z,t) 


§1 

§2 


§3 


—  u 

where  p  =  pA,  m  =  pAU,  e  =  pA(e  +  -j-)  are  conservation  variables  and 

P  =  PA,  h  =  e  +—  .  With  the  equation  of  state  for  perfect  gas  and 

P 

calorically  perfect  gas  assumption,  the  flux  vectors  f  can  be 
rewritten  in  terms  of  the  conservation  variables: 
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Hence  the  Jacobian  matrix 


r  o 


l^o2 

-y  H  -  (1-^u3 


i  o 

(3-Y)U  Y-l 

t^+T(1"y)ij2  yU 


The  complexity  of  this  Jacobian  matrix  makes  it  rather  tedious  to 
compute  the  eigenvalues.  It  is  easier  to  work  with  the  nonconserva¬ 
tive  or  primitive  variables  u  =  j^p  u  pjT.  We  can  rewrite  the 
equation  (A-l)  as 

Li  +  Ni  +  «-  °  <A-2> 

where 
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l _ 
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'  n  ' 

Multiplying  the  equation  (A-2)  by  the  inverse  matrix  L-"*",  we  convert 
it  into  a  nonconservative  form 


9u 
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(A-3) 
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U  p  0 

0  u  ~ 
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0  YP  U 


Si 

u 
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~  gl  + 

~s2 
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+  (1-y)  Ug2  +  (Y-l)g 

~h  ~ 
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-§3- 


and  yP  =  pc  where  c  =  ^yRT  is  the  speed  of  sound.  The  non¬ 
conservative  form  matrix  M  is  related  to  the  Jacobian  matrix  M 
by  the  similarity  transformation 


m  =  l_1ml 


Therefore  M  and  M  have  the  same  eigenvalues. 

It  is  easy  to  derive  the  eigenvalues  of  M,  which  are  U, 
U+c,  U-c.  The  corresponding  eigenvectors,  written  as  column 
vectors,  form  the  matrix 
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_  1  ~ 

such  that  T  M  T  =  A 
where 
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Multiplying  the  equation  (A-3)  by  T  ^  and  using  the  similarity 
relation  mentioned  above,  the  following  equation  is  obtained 


T_1  H  +  at-1  It  +  t_1  s  =  °  (a~4) 

This  equation  can  be  written  in  the  scalar  form  as 
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These  are  called  the  compatibility  equations  in  the  characteristic 
directions.  The  characteristic  directions  S^,  S2,  S3  are 
defined  by 

dq  =  ii  +  u  h  ■  then  if  -  D  alon«  si 

i  +  (U+c)  then  gf  -  U+c  along  S2 

dq  "  Tt  +  (U-o)  Tz •  then  if  "  U‘c  alonS  S3 


(A-6-1) 

(A-6-2) 

(A-6-3) 
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which  are  called  the  characteristic  equations.  Thus,  in  terms 
of  the  characteristic  coordinates,  the  compatibility  equations 
become 

c2  -  5^  +  c2  h  -  I3  -  0  alonS  S1  (A-7-1) 

(A-7-2) 
(A-7-3) 

It  may  be  noted  that  along  each  of  these  characteristic  curves, 
the  compatibility  equations  are  ordinary  differential  equations 
and  the  parameters  S-^,  S£,  can  be  considered  as  time  t 
itself.  Thus  we  can  rewrite  the  compatibility  equations  along 
the  characteristics  as  follows: 

(A-8-1) 
(A-8-2) 
(A-8-3) 
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APPENDIX  B 

BOUNDARY  CONDITIONS  BASED  ON  THE 
METHOD  OF  CHARACTERISTICS 


1.  Inlet  Boundary 

As  mentioned  in  the  main  text,  it  is  necessary  and  sufficient 
to  impose  two  boundary  conditions  at  the  inlet  boundary  for  sub¬ 
sonic  inflows.  In  the  present  case  the  total  pressure  P^  and  the 
total  temperature  T  are  prescribed  at  the  inlet.  For  determining 
all  the  variables  at  the  boundary  point,  however,  we  need  another 
equation  which  is  supplied  by  the  compatibility  equation  along 
the  characteristic  direction  . 


Q 

* 


\ 

\ 


At 

V 


\  Ay 

\aru-c 
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-x 


P  a  2 


J2‘ 


A 


In  general,  the  location  of  P  and  the  variables  at  Q  will 
have  to  be  solved  iteratively.  As  an  initial  iteration,  we 
approximate  the  characteristic  by  a  straight  line  with  the 
slope  evaluated  at  P.  (U-c)  at  P  is  approximated  by  interpola¬ 
tion  of  its  values  at  points  1  and  2.  Location  of  P  is  thus 
obtained  by  solving  the  following  three  linear  equations: 
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AL  +  A2  =  A  (B-l-1) 

Al  =  -  (U-c)p  At  (B-l-2) 

(U-c)p  =  [(U-c)1  A2  +  (U-c) 2  A1]/A  (B-l-3) 

Then  all  the  necessary  thermodynamic  variables  at  P  can  be 
determined  by  linear  interpolation  between  points  1  and  2. 

The  compatibility  equation  (A-8-3)  is  integrated  from  the  point 
P  to  the  point  Q  using  a  forward  difference  approximation,  which 
can  be  written  as 


F(Uq)  =  -  pc  (UQ  -  Up)  +  (PQ  -  Pp)  -  pc  At  g2+At  g3/A  =  0 

(B-2) 


where  Uq  and  Pq  are  unknowns  and  p  and  c  are  calculated  at  P  in 
the  initial  iteration.  Since  the  total  pressure,  Pt  and  the 
total  temperature,  Tt  are  prescribed  at  the  point  Q  we  can  express 
Pq  in  terms  Uq  by  the  equation 


u: 


Zc 


T 

P  t- 


Y 

Y-l 


(B-3) 


Using  this  equation,  we  can  eliminate  Pq  from  the  difference 
equation  (B-2)  and  obtain  the  following  nonlinear  equation  for  Uq. 


F(Uq)  =  -  pc (UQ  -  Up)  +  Pt 


-  Pp  -  pc  At  g2  +  At  g3/A  =  0 


(B-4) 


Differentiating  F,  we  can  write 
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The  equation  (B-4)  is  solved  by  the  Newton-Raphson  iteration: 


u(v+l)  =  u(v)_  F 

UQ  UQ  f 


(UQ(v))/F'(U<v)) 


(B-6) 


with  the  initial  guess  Uq°^  =  U-^ .  This  is  the  inner  iteration 
within  the  outer  iteration  started  earlier.  The  converged  value 
for  Uq  is  the  first  iterate  for  U  at  Q.  Further  iterations  are 
carried  out  to  satisfy  the  compatibility  equation  to  sufficient 
accuracy.  This  can  be  done  by  modifying  the  equation  (B-l-2) 
as 


Aq  =  -  Y  [  (U-c)p  +  (U-c)Q]  At  (B-7) 

where  the  latest  values  at  P  and  Q  are  used  in  each  iteration. 
Location  of  P  is  updated  in  each  iteration  and  in  the  equations 
(B-4)  and  (B-5)  p  and  c  are  calculated  by  the  average  of  their 
values  at  P  and  Q: 

P  =  y  (pp  +  Pq)  and  c  =  y  (°p  +  cq)  (B-8) 

This  iterative  process  is  continued  till  the  values  of  Uq 
converge.  This  convergence  is  usually  obtained  within  2  or  3 
iterations  at  each  time  step.  Once  Uq  is  known,  all  the  other 
variables  at  Q  can  be  computed. 

2.  Exit  Boundary 

For  subsonic  outflow,  only  one  boundary  condition  can  and 
should  be  imposed  at  the  exit  boundary.  The  static  pressure  P 
is  prescribed  at  the  exit  in  the  present  case.  The  other  vari¬ 
ables  at  the  exit  boundary  point  are  determined  by  using  the 
compatibility  equations  along  the  characteristic  directions 
and  S2 . 
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a.  Integration  along  the  character ic  direction  S^: 


In  a  manner  similar  to  the  procedure  used  at  the  inlet 
boundary,  the  following  three  linear  equations 


A1 

+  A2  -  A 

(B-9-1) 

A2 

=  UgAt 

(B-9-2) 

US 

=  [UI_1  A2  +  Ux  A1]/A 

(B-9-3) 

are  solved  to  determine  the  location  of  S  and  the  thermodynamic 
properties  of  the  flow  at  the  point  S  are  determined  by  a  linear 
interpolation  between  the  points  I  and  (1-1) .  The  compatibility 
equation  (A-8-1)  is  integrated  by  a  forward  difference  approxi¬ 
mation  and  the  density  at  the  boundary  point  0  is  calculated  by 

P0  =  Ps  +  [PQ  -  Pg  -  At  (c2  gx  -  g3)/A]  /c2  (B-10) 

where  Pq  is  known  from  the  imposed  boundary  condition  and  all 
the  coefficients  are  calculated  at  the  point  S. 
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b.  Integration  along  the  characteristic  direction  S^: 


The  same  procedures  as  in  case  (a)  is  followed  except  that  U  in 
the  equations  (B-9-2)  and  (B-9-3)  is  replaced  by  U+c .  The  com¬ 
patibility  equation  (A-8-2)  is  integrated  to  obtain  U  at  the 
point  0  as 

u0  =  UR  -  l>0  -  PR  +  At  (pc  g2  +  g3/A)]/(pc)  (B-ll) 

The  entire  process  of  (a)  and  (b)  can  now  be  repeated  by  replacing 
U  in  the  equation  (B-9-2)  and  c  in  the  equation  (B-10)  and  pc  in 
the  equation  (B-ll)  by  their  respective  averages  at  the  points 
0  and  S  or  0  and  R.  Convergence  is  usually  achieved  in  2  or  3 
iterations.  All  the  variables  necessary  at  the  exit  boundary 
point  0  can  now  be  determined. 
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c.  Station  locations 


d.  Elemental  control  volume 
Figure  1.  Concluded. 
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Figure  4.  Model  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  87  percent  corrected  speed  (Runge-Kutta  Scheme). 


AEDC-TR-82- 1  6 


PRESSURE  RATIO  ACROSS  THE  COMPRESSOR 


5.  70, 


5.  60 

ce  5.  50 

o 

GO 

£3  5-  40 

Cl 

§  5.30 

C_J 

uj  5.  20 

i  1 

;  0^Exper 

i mental  Surge 

Point 

\ 

i — 

oo  5.  10 

GO 

S  5.  m 

<c 

2  4.  90 

1 - 

<C 

^  4.  80 

\ 

> 

\ 

1 

LU 

Cel 

=  4,  70 

GO 

OO 

LU 

gf  4.  60 

4.  50 

I 

4.  40 

s 

Q 

(S3 

C3 

51 

53 

51 

53 

Q 

Q 

(3 

53 

(3 

(3 

• 

• 

a 

0 

a 

a 

0 

CM 

cn 

in 

CO 

t". 

00 

CD 

CD 

CD 

CO 

CD 

CO 

CO 

PERCENT  CORRECTED  COMPRESSOR  INLET  AIRFLOW 


Figure  5.  Model  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  87  percent  corrected  speed  (JRS  Scheme). 
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Figure  6.  Model  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  87  percent  corrected  speed  (MacCormack  Scheme). 
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Figure  7.  Variation  of  the  compressor  flow  parameters  with  time  till  surge  point  -  87  percent 
corrected  speed  (Runge-Kutta  Scheme). 
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Figure  8.  Variation  of  the  compressor  flow  parameters  with  time  till  surge  point  -  87  percent 
corrected  speed  (JRS  Scheme). 
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Figure  9.  Variation  of  the  compressor  flow  parameters  with  time  til  surge  point  -  87  percent 
corrected  speed  (MacCormack  Scheme). 
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Figure  10.  Individual  stage  pressure  variai 
speed  (Runge-Kutta  Scheme). 
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Figure  13.  Individual  stage  airflow  variations  just  before  surge  -  87  percent  corrected 
speed  (Runge-Kutta  Scheme). 
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Figure  14.  Individual  stage  airflow  variations  just  before  surge  -  87  percent  corrected 
speed  (IMS  Scheme). 
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Figure  15.  Individual  stage  airflow  variations 
speed  (MacCormack  Scheme). 
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Figure  16.  Model  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  102  percent  corrected  speed  (Runge-Kutta  Scheme). 
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Figure  17.  Mode!  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  102  percent  corrected  speed  (JRS  Scheme). 
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Figure  18.  Model  prediction  of  the  compressor  performance  under  loading  till  surge  point 
at  102  percent  corrected  speed  (MacCormack  Scheme). 
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Figure  20.  Variation  of  the  compressor  flow  parameters  with  time  till  surge  point  -  102 
percent  corrected  speed  (IMS  Scheme). 
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Figure  21.  Variation  of  the  compressor  flow  parameters  with  time  till  surge  point  -  102 
percent  corrected  speed  (MacCormack  Scheme). 
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e  22.  Individual  stage  pressure  variati 
speed  (Runge-Kutta  Scheme). 
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Figure  23.  Individual  stage  pressure  variations  just  before  surge  -  102  percent  corrected 
speed  (JES  Scheme). 
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Figure  24.  Individual  stage  pressure  variations  just  before  surge  -  102  percent  corrected 
speed  (MacCormack  Scheme). 
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